Spectral methods for the wave equation in second-order form 
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Current spectral simulations of Einstein's equations require writing the equations in first-order 
form, potentially introducing instabilities and inefficiencies. We present a new penalty method for 
pseudo-spectral evolutions of second order in space wave equations. The penalties are constructed 
as functions of Legendre polynomials and are added to the equations of motion everywhere, not 
only on the boundaries. Using energy methods, we prove semi-discrete stability of the new method 
for the scalar wave equation in flat space and show how it can be applied to the scalar wave on a 
curved background. Numerical results demonstrating stability and convergence for multi-domain 
second-order scalar wave evolutions are also presented. This work provides a foundation for treating 
Einstein's equations directly in second-order form by spectral methods. 
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I. INTRODUCTION 

Recent advances [1-3] in numerical simulations of black 
holes in general relativity have led to many interesting 
results. Most of these simulations have been carried 
out with finite-difference methods. However, the vac- 
uum Einstein equations have mathematically smooth so- 
lutions (unless pathological coordinates are chosen) . Ac- 
cordingly, one expects that spectral methods should be 
optimal in terms of efficiency and accuracy. 

Einstein's equations are a hyperbolic system involv- 
ing second derivatives in space and time. However, 
the numerical solution of hyperbolic systems using spec- 
tral methods is normally performed with a fully first- 
order formulation, even when the equations are natu- 
rally higher order. Reducing the order of the equations 
is usually achieved by introducing new variables defined 
as first-order time or space derivatives. The basic im- 
petus for this first-order reduction is that there exists 
a well-established body of mathematical literature for 
first-order hyperbolic systems [4-6] , which includes meth- 
ods for analyzing the well-posedness of the equations and 
the proper way to impose stable boundary conditions in 
terms of characteristic variables. 

The obvious disadvantage of the first-order reduction is 
the introduction of additional variables, whose definitions 
(at least for spatial derivatives) become constraints the 
solution must satisfy and thus new possible sources of 
instability in the system. Furthermore, each new variable 
must be evolved, increasing the number of equations and 
the computational cost of the simulations. In some cases, 
this can be a substantial increase. 

Successful simulations of Einstein's equations using 
spectral methods have thus far been implemented only 
as first-order reductions of the second-order system [7, 8] . 
In the case of the generalized harmonic form of the equa- 
tions, the reduction to first order in space proceeds by 
introducing 30 additional variables, more than doubling 
the number of equations and constraints in the system [7] . 
These simulations typically require significant computa- 



tional time, upwards of a hundred CPU-weeks for high 
resolution runs [8] . 

A first order in time, second order in space system 
has the potential to reduce the constraint-violating in- 
stabilities and the computational expense of the simu- 
lations. However, the mathematical knowledge under- 
lying the proper formulation for such systems is much 
less developed. Recently, Gundlach and Martin-Garcia 
have proposed and analyzed definitions of symmetric hy- 
perbolicity for a general class of second order in space 
systems [9, 10]. They have also shown how one may de- 
fine characteristic modes in the second-order system and 
thereby formulate stable boundary conditions at the con- 
tinuum level. 

There still remains the problem of how to impose the 
boundary conditions in the discrete system (using spec- 
tral methods). Even for the simplest representative hy- 
perbolic system, the second order in space wave equa- 
tion, naive attempts to impose boundary conditions in 
the same way as in a first-order formulation generally 
fail. The difficulty is not due solely to the presence of sec- 
ond derivatives. For example, methods exist for treating 
the second order spatial derivatives in the Navier-Stokes 
equations directly using spectral methods [11-13]. How- 
ever, these techniques do not apply to the wave equation, 
as the characteristic structure is fundamentally differ- 
ent. In this work, we present a new method for imposing 
boundary conditions in the second-order wave equation 
that is robust, stable, and convergent. 

Since the generalized harmonic form of Einstein's equa- 
tions appears as ten nonlinear coupled wave equations, 
this work provides a foundation for solving Einstein's 
equations directly in second-order form using spectral 
methods. This application will appear in a subsequent 
paper [14]. It is likely that the work presented here 
will also allow other formulations of Einstein's equa- 
tions, such as the BSSN (Baumgarte-Shapiro-Shibata- 
Nakamura) formulation, to be treated by spectral meth- 
ods without reduction to first-order form. 

In Section II A we review a typical spectral method for 



evolving the fully first-order form of the one-dimensional 
wave equation. We review how boundary conditions can 
be imposed using penalty methods [15] and how stability 
of the system can be analyzed with energy methods [4, 5]. 
In Section II C we present the new penalty method for the 
one-dimensional second order in space wave equation and 
prove stability of the system using energy arguments. In 
Section III we generalize the method to three dimensions, 
and in Section IV we apply the method to the case of a 
scalar wave on a curved background. 



II. ONE-DIMENSIONAL WAVE EQUATION 

We begin with the one-dimensional wave equation 

i>^r, (1) 

where V' = i^ix, t). Here, dots denote differentiation with 
respect to i, while primes denote differentiation with 
respect to x. We will first review a typical first-order 
pseudo-spectral method for evolving this equation before 
discussing the second-order formulation. 

A. First-Order System 

The wave equation in Eq. (1) reduces to first order by 
introducing the variables -k and (^, where 






(2) 
(3) 



The negative sign in the first equation is purely a matter 
of convention. The first-order system is thus 



■k — 






(4) 
(5) 
(6) 



Equation (4) is just the definition of tt, while the defi- 
nition of (/) in Eq. (3) amounts to the addition of a con- 
straint C = to the system, where 



C = ^' 



(7) 



The characteristic variables U and speeds A for this sys- 
tem (see e.g. Ref. [4]) are 






A = 0, 
A = ±l. 



(8) 
(9) 



Here, n^ is the unit outgoing normal vector to the bound- 
ary, which in one dimension is just n^ = ±1. With this 
definition, U- is incoming (A < 0) at each boundary. 

For a symmetric hyperbolic system on a domain D, with 
boundary dfl, there exists a (not necessarily unique) con- 
served, positive definite energy 



E= / edV, 
In 



(10) 



which is conserved in the sense that 
e = d,F\ 



(11) 



Accordingly, the time derivative of the energy is given by 
the fiux through the boundary. 



E 



F^dA, 



on 



(12) 



where F" = n^F*. Note that for general quasi-linear 
systems, the energy is only strictly conserved when co- 
efficients in the equations are approximated as constant 
and lower-order terms are neglected. 

For the one-dimensional wave equation in Eqs. (4)-(6), 
the energy density is 



1 



(^2+02). 



Using Eqs. (11), (5), (6), and (9), we get 



n 



F"" = -TT (t) = '—{Ul -Ul). 



(13) 



(14) 



If we consider our domain to be the interval [—1, 1], then 
E-\T.(U'--Ul). (15) 

x = ±l 

For well-posedness and stability, one requires that the 
growth of the energy be controlled by specifying bound- 
ary conditions for the positive terms in E. Therefore, 
a boundary condition must be supplied on the incoming 
mode U-. For example, with a homogeneous condition 
specifying C/_ — (or more generally U- = k t/+ for 
|k| < 1), it follows that E < 0. Together with the positive 
definiteness of the energy, this ensures that the system 
is stable. If instead the incoming mode is a prescribed 
function [/_ = /, then wc still have stability in the sense 
of controlling the energy with a bound that involves /. 

The definition of energy given by Eq. (13) is not 
unique, but was motivated in part by a desire to obtain a 
sharp energy bound. For example, we could have defined 
the energy density with a term a^V^ ^^ 



e=l(a'^'+n' + c^'). 



(16) 



In this case, we would obtain the additional term in E 



a'^xPndV <- / (aV^ + tt^) dV, (17) 



where the inequality follows from the relation 2uv < u^ + 
v^ for any (real) u, v. We would thus arrive at the weaker 
estimate 



^^iE(^--f^+ 



aE. 



(18) 



x=±l 



With a condition on the incoming mode C/_ to control 
the boundary term on the right-hand side, the system 
is stih weh-poscd in this case [4], but we are unable to 
prove stability in the sense that _E < 0. 

In the semi-discrete problem, one considers the dis- 
cretization in space but not time. An effective way to im- 
pose boundary conditions in the semi-discrete system is 
to add appropriate penalty terms to the equations on the 
boundaries. Such a penalty method thus imposes con- 
ditions "weakly" — that is, approximately, without com- 
pletely replacing the equation of motion on the bound- 
ary [15]. Heuristically, the rationale for this is that it is 
not necessary to enforce exact boundary conditions on 
approximate solutions. All that is required is for the 
discrete solution to converge to the continuum solution 
with the correct boundary conditions as the resolution 
is increased. We find that these methods generally yield 
superior accuracy and convergence while also providing 
a simple way to impose arbitrary boundary conditions. 

The penalties are added to the equations on the bound- 
ary in the form {U^'^ — U-), so that if the boundary 
condition [/_ = C/^"-' is satisfied then the penalties van- 
ish. The appropriate penalty (up to an overall coefficient) 
for each equation can be found by projecting the bound- 
ary conditions in terms of characteristic variables to fun- 
damental variables [16]. In other words, we first trans- 
form to characteristic variables in the first-order system 
of Eqs. (4)- (6) on the boundary and add penalties: 

1 



u^^--{u+ + u^), 




(19) 


U+ = -n-U[, 




(20) 


U^ = +n^U'_+c{U^^ - 


-U-), 


(21) 



Suitable values for the penalty factor c in Eqs. (23)- 
(24) can be determined from a semi-discrete energy anal- 
ysis, which we will now show. For case in obtaining an- 
alytical results, we choose Gauss- Legendre-Lohatto collo- 
cation points (see Appendix A for details). The basis 
functions for this choice are the Legendre polynomials 
Pn{x) on [—1, 1]. We begin by writing the semi-discrete 
energy corresponding to Eq. (13): 



E 



[(7r,^) + (<^,</))] 



(25) 



where 



represents a discrete inner product, as in 



N 



(ti", tt) = ^ Wj 7r| 



(26) 



Here tt^ are the grid values of the function tt, and uji are 
the quadrature weights (see Appendix A). Taking the 
time derivative of the semi-discrete energy in Eq. (25), 
we obtain 



I A^ C 



(27) 



where we have used summation by parts (the discrete 
analogue of integration by parts) in the first term and 
introduced the notation SU^ = C/^^ — U- in the penalty 
terms. The first term in Eq. (27) is similar to the con- 
tinuum result: 



-^^^l"o = iE(^--f^+)- 



(28) 



=0.N 



where c is an undetermined constant. Only the equa- 
tion for U- has a penalty term, since there is no bound- 
ary condition on U^ or C/_|_. We then transform back to 
fundamental variables to obtain the first-order equations 
with penalties: 

V'« = - ^«, (22) 

7T. = - 0^ + I (<5,o + <5.N)(t/^^ -U-), (23) 

<i>^^-<-^ n^S.o + SmW^"" - U-). (24) 



Here we have explicitly denoted grid values with a sub- 
script i. For a pseudo-spectral method, one chooses the 
nodes of a Gaussian quadrature rule as collocation points. 
The A^ + 1 grid points Xi run from xo = — 1 to xjv = +1. 
Differentiation is implemented by matrix multiplication, 
as in tt' = ^1 ■ -D,-, tt, , where D)-' is the first-order dif- 
ferentiation matrix. The Kroncckcr delta terms dio + SiN 
indicate that the penalties are applied only on the bound- 
aries at i = 0,N. The penalty coefficients should satisfy 
c — > oo as A^ — i> cxD, in order to ensure that the continuum 
equations and boundary conditions are recovered in this 
limit [15]. 



Evaluating the discrete inner products in the last two 
terms in Eq. (27) yields 



-C^pcnaltii 



f T. ''-«-. 



(29) 



i=O.N 



where we have written uj for the quadrature weight ujq 
lun at X = ±1. Noting that 



1 



_(^BC2 



U-5U- 
we put things together to find 



t/2 _ SU^) 



(30) 



^=1 E [{^-cu;)Ul-Ul + cu;iU^^^-SUl)]. (31) 



i=0,N 



The condition on the penalty factor c for stability de- 
pends on the boundary condition we impose on U_. 
Requiring E < 0, we find: 



U: 



BC 







1 

=> c > — , 

UJ 



U^^ = kU, 



(jj K^ 



> C > -, 

UJ 



(32) 
(33) 



where |k| < 1. The strictest condition is obtained by 
insisting that the energy be bounded by the continuum 
energy in Eq. (15) for arbitrary U^'~^: 



E < E, 



continuum 



_ 1 



(34) 



The situation is sHghtly different when considering the 
semi-discrete energy for a multi-domain problem. For 
example, suppose we consider the interval [—2,2] with 
an inner boundary at a; = 0. The energy calculation 
up to Eq. (31) is identical on each subdomain. The key 
difference is that now the incoming mode at the interface 
boundary is supplied by the adjacent subdomain. If we 
denote the intervals [—2, 0] and [0, 2] with subscripts 1 
and 2, respectively, then at a; — 0: 



C/f- = U2+, 



c/P? 



Uu 



The terms in _E at a; = are then 

r2 /I „ , .\7r2 



(1 - cuj)U^_ - (1 - cw)t/i% - cu}5U^_ 
+ {1- cuj)Ul_ - (1 - ccj)t/|+ - cu}SUi_ 



(35) 
(36) 



(37) 



This quadratic form is negative semi-definite if and only 
ii c = 1/lu. On a multi-domain problem, the value of c 
at an internal boundary required for stability is there- 
fore fixed, regardless of the boundary condition imposed 
at the external boundaries. This analysis assumes that 
the penalties at the interface boundary enforce conditions 
only on the incoming modes. It is also possible to penal- 
ize arbitrary combinations of the variables at interfaces 
and thereby to obtain different stability conditions (see 
e.g. Ref. [17]), but we do not consider this refinement 
here. 

On an arbitrary domain, the definition of the discrete 
inner product is modified. For instance, if we consider 
a one-dimensional domain il with a coordinate mapping 
fx : [—1,1] -^ ^, then the Jacobian of the mapping is 
inherited from the continuum inner product: 



if, 9) 



N 

E 



^ifigifj''f 



(38) 



points. In particular, a penalty that is applied only on the 
boundary of a Legendre grid as in Eqs. (23)-(24) would 
in general be non-zero everywhere on a Chebyshev grid. 
In practice, the system works well even without modifica- 
tion on a Chebyshev grid by simply using Eqs. (22)-(24) 
as derived for a Legendre grid and letting the index i 
represent the Chebyshev grid points. Chebyshev stabil- 
ity of penalty methods is proved in Ref. [15] for simple 
cases, and proofs of Chebyshev stability for more general 
hyperbolic problems are reviewed in Ref. [19]. 

It is also worth noting that stability conditions derived 
from strict energy arguments can generally be relaxed to 
a degree. The penalty factor c, which was found to be 
l/cj = N{N + l)/2 for Legendre methods, can be opti- 
mized by trial and error to maximize efficiency and ob- 
tain the least restrictive Courant-Friedrichs-Lewy (CFL) 
condition while maintaining stability. This is discussed, 
for example, in Ref. [11]. 

Stability of the fully discrete problem can be explored 
by examining eigenvalues. The entire system is written 
in a form suitable for passing to an explicit time-stepping 
algorithm, as in y = Ay, where the vector y represents 
the grid values of all the fields. The eigenvalues of the 
matrix A can then be plotted in the complex plane and 
compared with the stability region of the time-stepper. 
In general, positive real parts of the eigenvalues imply in- 
stability, while the spectral radius (maximum amplitude 
of eigenvalues) is inversely proportional to the maximum 
allowed time-step [4] (the exact relation depends on the 
time-stepping algorithm being used). 

A typical eigenspectrum for the system in Eqs. (22)- 
(24) on two subdomains is shown in Fig. 1. Curiously, the 
large amplitude conjugate pair of eigenvalues on Cheby- 
shev points is absent on the Legendre grid. This implies 
that there is a less restrictive CFL condition for the sys- 
tem on Legendre grid points, and this is indeed the case 
for this particular system. As is discussed elsewhere, it is 
unlikely that this difference carries over to more general 
systems [5]. For instance, we find no significant differ- 
ence in time-stepping conditions on Chebyshev or Leg- 
endre grids for the three-dimensional wave equation (in 
fiat or curved space). It is also worth noting that eigen- 
value stability is insufficient to prove that the system is 
actually stable and convergent, but it is suggestive [5]. 



Since the penalty terms in Eq. (27) contain Kronecker 
deltas that pick out specific terms from the sums, the 
values for c we obtain would need to be modified by a 
Jacobian factor: c -^ c/ fx' . For simplicity, we will assume 
that the domain is the fundamental interval [—1, 1] unless 
otherwise stated, so that no Jacobians are needed. 

Although we performed the semi-discrete energy anal- 
ysis on Legendre grid points, this is not a limitation. 
One could implement the system on Gauss-Chebyshev- 
Lobatto points using, for example, the Chebyshev- 
Legendre method [18]. With this method, the equations 
are implemented on a Chebyshev grid by interpolating 
the (Legendre-grid) penalty functions to the Chebyshev 



B. Second Order in Space 

The first order in time, second order in space formulation 
of the one-dimensional wave equation in Eq. (1) is 

V- = - TT, (39) 

7T = - ?A". (40) 

The characteristic variables are the same as those of the 
first-order reduction with the replacement (/)—>■ ^': 

C/^=^, A = 0, (41) 

C/±=7r±n^V', A = ±l, (42) 
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FIG. 1: Left: Eigenvalues in the complex plane of the first-order system, Eqs. (22)-(24). 
unstable second-order system, Eqs. (45)-(46). Both plots are for a two-domain problem on 
a; = 0, penalty factors c — N{N + l)/2, outer boundary conditions f/j? — 0, and A^ -I- 1 
Results for Legendre- and Chebyshev-Lobatto grids are shown for comparison. 



Right: Eigenvalues of a typical 
-1,1], with an inner boundary at 
= 11 grid points per subdomain. 



The energy and flux are the same as well: 
1 



e=-{n^+i^'') 



^=1 E(f^--^+)- (43) 



x=±l 



The difhculty arises in the semi-discrete second-order sys- 
tem when trying to find appropriate penalties by project- 
ing from characteristic variables, as was done in the first- 
order case. The boundary condition SU- = C/^'^ — C/_ = 
is now a differential as opposed to an algebraic condi- 
tion: 



n^V' = U- 



BC 



(44) 



One therefore obtains a condition on tp' at the boundary, 
but not on tp itself (there is no boundary condition on 
U^). The system one arrives at by naively following the 
same procedure as in the first-order case is 



ipi 



i^i +c{Sio + SiN)SU-. 



We might try applying a penalty to Eq. (45) also: 

TTj = - V'-' + C2 ((5jO + 5in) SU-. 



(45) 
(46) 



(47) 
(48) 



These equations are generally unstable, particularly 
when evolved on multiple subdomains with at least one 
interface boundary. The error in ip tends to grow ex- 
ponentially, ruining the evolutions within a few hundred 
crossing times. The penalty factors ci,C2 can be fine- 
tuned by trial and error to obtain approximately sta- 
ble evolutions in some cases, but not robustly so. Fig- 
ure 1 shows a typical eigenspectrum for the system in 
Eqs. (45)-(46) on two subdomains. The eigenvalues with 
positive real parts clearly indicate instability. 



C. Second-Order Penalty Method 

We will now derive a way to impose penalty boundary 
conditions in the second-order system that yields a ro- 
bust, stable result. For the semi-discrete problem, we 
once again choose Gauss-Legendre-Lobatto collocation 
points. We begin by writing the second-order equations 
in the form 



7T, = - V" + q, 



(49) 
(50) 



where p and q represent as yet undetermined penalties. 
The semi-discrete energy is 



1 



Taking the time derivative, we find 



(51) 



E - -V'-Tr.f + V-ftf + (7r,q) - (V^",p), (52) 

where we have used summation by parts in the first two 
terms. The first term is like the continuum result: 






lo 



(53) 



Since the projection of boundary conditions from charac- 
teristic to fundamental variables is unambiguous in the 
variable tt, we will write the penalty q as in Eq. (46): 



— (Sio 

UJ 



iN 



5U-, 



(54) 



where a is an undetermined constant and w is the quadra- 
ture weight at a; = ±1. The factor l/u is explicitly writ- 
ten in anticipation of its cancellation when evaluating the 
third term in Eq. (52): 



(tt, q) — aTiQ (5C/° + aiTM5U_ 



N 



(55) 



If we also choose the penalty p in Eq. (49) to have a 
similar value on the boundary 



p — a SU- , 



(56) 



then the second and third terms in Eq. (52) combine to 
form the expression 



1=0,^ 



BC2 



Ul -6Ul). (57) 



Note that we do not define p on the boundary with a 
factor of l/oJ, because the second term in Eq. (52) arises 
out of summation by parts as opposed to being picked 
out from the discrete sum by a Kronecker delta. 

The stumbling block in this energy analysis is the last 
term in Eq. (52), whose appearance is inevitable because 
of the derivatives in the definition of energy in Eq. (51). 
Such a term did not arise in the first-order energy esti- 
mate, precisely because the first-order energy in Eq. (25) 
did not contain any derivatives. Fortunately, it turns out 
we can eliminate the inner product {'ip",p) by allowing 
the penalty p to be non-zero throughout the domain and 
by constructing it to be orthogonal to ip" . 

The scalar field "0 in the semi-discrete solution is an 
interpolating iV*'^-order polynomial [6]. Therefore, ip" is 
an iV — 2 order polynomial, and the product ip"p is at 
most a polynomial of order 2N — 2. It follows that the 
quadrature integral is exact: 
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FIG. 2: Eigenvalues in the complex plane of the second-order 
system in Eqs. (64)-(65) on two subdomains covering [—1, 1], 
with an inner boundary at a; = 0, outer boundary conditions 
U- = 0, and A'^ -I- 1 = 11 grid points per subdomain. Re- 
sults for Legendre- and Chebyshev-Lobatto grids are shown 
for comparison. 



which is just like Eq. (31) for the first-order system with 
2a o cw. The conclusions reached previously for c there- 
fore carry over: a multi-domain problem with arbitrary 
outer boundary conditions is stable only if a = 1/2. The 
second-order system with penalties is thus 



(V'",p)= r{x)p{x)dx. 



(58) 



This inner product will automatically vanish if the 
penalty p is a linear combination of the Legendre poly- 
nomials Pn{x) and PN-iix), which are orthogonal to 
any polynomial of degree A^ — 2 or less. We are therefore 
provided with two degrees of freedom for constructing 
the function p, which is sufficient to satisfy the boundary 
values defined in Eq. (56). We make use of the follow- 
ing polynomials, constructed to take the values 0, 1 at 
a; = ±l: 



1 



fix)^-{-ir[PN{x)-PN-l{x)], 

9{x) ^-[PNix)+PN-i{x)]. 
If we now define the penalty p to be 

P^Pafix) +PN9{x), 



(59) 
(60) 

(61) 



where po and pn represent the endpoint values of 
Eq. (56), then the penalty function p will have the correct 
boundary values while also satisfying 

(Vy',p)-0. (62) 

Putting things together, we now obtain 

^ = i E [i^-2a)Ul-Ul + 2a{U^''^-SUl)], (63) 

i=0,N 



^, = - TT, - - [fix) SUl + g{x) SU"!] , (64) 
ir^=-^'l + ^[S^o5Ul+5,N5U'!]. (65) 



Of course, one needs to be concerned not only with stabil- 
ity^ but also consistency — that is, the system should re- 
produce the continuum equations in the limit as A^ — ?> oo. 
The penalty p on the ijj equation in Eq. (64) is applied 
throughout the domain and not only on the boundaries. 
Moreover, it does not scale as iV^, so consistency might 
seem dubious. However, the penalty on -k in Eq. (65) 
does scale as N"^ and is applied only on the boundaries. 
Therefore, the condition 5U- -^ on the boundary as 
iV — ^ oo is enforced. This also implies p — ^ in turn, so 
consistency follows. 

Although the second-order energy argument was per- 
formed on Legendre points, the equations can be imple- 
mented on any grid, just as in the first-order system dis- 
cussed in Section II A. Eigenvalues of the fully discrete 
system imply stability here as well, as shown in Fig. 2 for 
a representative two-domain problem. The spectral ra- 
dius is somewhat larger than in the first-order spectrum 
in Fig. 1. However, we find that differences in CFL con- 
ditions essentially disappear for more general systems, 
including the three-dimensional wave equation in flat or 
curved space. 



III. THREE-DIMENSIONAL WAVE EQUATION 

Let us now consider the three-dimensional second order 
in space wave equation 



V" = -TT, 

TT = - did''ip. 



(66) 

(67) 



The characteristic modes and speeds of this system are 



U^ =ip, A = 0, 

L/± =7r±n'9,V, A = ±l, 

C/f = diip — tiin^djip, A = 0, 



(68) 
(69) 
(70) 



where rf is the outward-directed unit normal to the 
boundary. These are the same as the characteristic 
variables of the first-order system obtained by defining 
(j>i = ditp. Usually, one thinks of characteristic variables 
as being defined only for first-order systems, but they 
can be generalized to second-order systems. One way 
to do this is to define the second-order modes as those 
combinations U of variables (tt, diip) that satisfy 



U = -Xn'd,U 



(71) 



where the dots represent derivatives transverse to n^ 
plus lower order terms [10]. As a consequence of this 
definition, the transverse derivatives Uf are automati- 
cally zero-speed modes (in fact, they can be given arbi- 
trary speeds). Moreover, the characteristic variables in 
Eq. (69) are unique only up to addition of these zero- 
speed modes. For example, we could redefine U± as 
U± + X'C/° for arbitrary (fixed) X*. As discussed in 
Ref. [10], this ambiguity is removed for a symmetric 
hyperbolic system by requiring the existence of a con- 
served energy that is quadratic in the modes. Here, that 
amounts to taking the definitions in Eqs. (68)-(70) as 
they are. The conserved energy density is 



e = i(^2^ava,7/.). 



(72) 



Note that this energy is indeed quadratic in terms of the 
characteristic modes: 



6=i(t/^ + C/^)+t/0»(7°. 



(73) 



In analogy with the one-dimensional case in Eq. (15), the 
flux is 



4 



dn 



{Ul - Ul) (fx, 



(74) 



where 9f7 represents the boundary of the domain. 

Now consider the semi-discrete problem in three- 
dimensions. We encounter a few issues in generalizing 
from the one-dimensional case. For one, if the bound- 
ary of the domain contains edges or corners, the normal 
vectors there (and hence characteristic modes) are not 



well-defined. For reasons that will become clear below, 
we resolve this ambiguity by defining the normal vectors 
as follows. We will use upper case N and lower case 
n to denote the unnormalized and unit normal vectors, 
respectively. For simplicity, suppose the domain il is a 
cube with x,y, z G [—1,1]. On boundary faces (codimcn- 
sion 1), one coordinate is fixed (e.g. the x = +1 face). 
We define face normals on a boundary with a fixed i*^ 
coordinate as 



N = LUj ujk n, 



(75) 



where ujj and ujk are the quadrature weights (see Ap- 
pendix A) corresponding to the two free dimensions, and 
n is the usual (Cartesian) unit normal vector in the z*^ 
direction. On edges and corners, the normal vector is 
defined to be the sum of the normals to the adjacent 
boundary faces. For example, the normal vector at the 
corner (x, y, z) = (1, 1, 1) is defined to be 

N ^ UJyUJzX + UJxUJzy + i^x^y2.. (76) 

The second-order system with penalty functions p, q is 



tp = - -K + p, 

n = — did'tp + q, 



(77) 
(78) 



where for conciseness we have suppressed indices repre- 
senting grid values (e.g. ip = ipijk)- Consider the semi- 
discrete energy 



s==i[(7r,^) + (av,ar 



(79) 



The computation of E proceeds analogously to the one- 
dimensional case, except for complications due to the cor- 
ners and edges. To see this, consider the following term 
that arises in taking the time derivative of Eq. (79): 



(aV,9i7r) = '^OJiLJjOJkd'^tpdlTT. 



(80) 



We use summation by parts in this expression and obtain 
three boundary terms — one for each I. For example, from 
I ~ z we get 



Each point on an edge receives a contribution from two 
such boundary terms, while points on corners get a con- 
tribution from all three. On the cube at the corner point 
(1, 1, 1), for example, the value obtained is 



{ui^Uly dzTp + L^x^z dyip + U!yU!z dxlp) TT. 



(82) 



We would like to be able to write this in terms of char- 
acteristic modes as 



N'ditp n 



INI 



{Ul - Ul) 



(83) 



and this is precisely the reason for the definition of nor- 
mal vectors on edges and corners given above. Thus, 
Eq. (80) can be written 

(aV,a.7r) ^^J2\N\{Ul - Ul) - {d^d.^,7:), (84) 
an 

where the sum is over all boundary points, including 
edges and corners. The magnitude of the normal vector 
N| encodes the appropriate quadrature weight factors 
for boundaries of any codimension. 

In a similar way, the terms in E due to the penalty p 
in Eq. (77) can be written 

(av,5.p) - Y. \^yd,^pp^{d'd^^p,p). (85) 

The penalty q in Eq. (78) is applied only on the boundary, 
where it takes the value 



1 INI 



dn 2 bJxLOyUJz 



-SU. 



(86) 



On a boundary face with fixed i"^-coordinatc, this re- 
duces to 



2uji 



(87) 



just as in the one-dimensional system. Assuming the 
boundary values of p satisfy 



P 



dn 



-Iju^ 



the penalty contributions to E combine to give 

penalties = ^J2 \^\U-SU- - {d'd^^P,p}. 



dn 



With everything included, the energy flux is 

^ = I E |N| (u^^^ -ul- 5Ul) - {d'd,^,p). 



(88) 



(89) 



(90) 



an 



We would like to eliminate the last term with an appro- 
priate choice of bulk penalty function p, as was done in 
the one-dimensional case. The most obvious generaliza- 
tion of the one-dimensional approach would be to con- 
struct p out of polynomials 6'„ satisfying {d'^di'ip, On) = 0. 
Unfortunately, this cannot be done. There are in general 
only about 2N'^ such functions 6'„ — not enough to sat- 
isfy (SN'^-\-2 boundary conditions (a proof is provided in 
Appendix B). 

Alternatively, one could seek a solution by allowing the 
penalty p to depend explicitly on the scalar field i/". One 
way of doing this is to split the offending inner product 
term of Eq. (90) into contributions from the boundary 
and the interior of the domain: 



(A^,p) = (A^,p) 



dn 



(Ai/.,p)| 



(91) 



where A-^ = d^diip. Considering the values of p on the 
boundary to be specified by Eq. (88), the first term on 
the right-hand side of Eq. (91) is fixed. We can then 
define p on the interior of the domain to be 



Pintc 



{Aip,p) 



an 



(AV',AV') 



AV', 



(92) 



provided A-0intcrior / 0. With this definition, the discrete 
sum over the interior cancels the sum over the boundary 
in Eq. (91), and the inner product (A^,p) vanishes. If 
Alp = 0, we cannot use Eq. (92), but in this case there 
would be no need since then (Ai/),p) = 0. One way this 
recipe could fail is if Aip vanishes on the interior of the 
domain but not on the boundary, although we believe this 
to be very unlikely in an actual numerical simulation. We 
find that while this method yields stability, it suffers from 
lack of convergence as resolution is increased. However, 
as we have not experimented extensively with approaches 
like this, further investigation may be worthwhile. 

Abandoning any explicit dependence on ip in the penal- 
ties, we seek instead to construct p so as to minimize the 
inner product {d^diip,p) of Eq. (90). It turns out this 
can be done by using a penalty constructed out of the 
same functions f,g defined in Eqs. (59)-(60) for the one- 
dimensional problem. Here we will give a summary of 
the result; a derivation is provided in Appendix C. We 
define one-dimensional functions /, g along each dimen- 
sion and write their grid values as fi — ,f{xi), fj — ,f{yj), 
fk = f{zk) (and similarly for g). Assuming the values of 
the penalty function p on the domain boundary dfl are 
given, the grid values on the interior of the domain are 



Pijk — POjk fi + PNjk gi + PiOk fj + 
— POOk fifj — PNOk gifj — . . . 

+ Pooo fifjfk + Pma gifjfk - 



(faces) 

(edges) (93) 
(corners) 



The bulk penalty picks up a contribution from each 
boundary face, edge, and corner. The assumption of a 
cubic domain is not a limitation, as it is straightforward 
to generalize this procedure to other domains. With this 
choice of penalty, E is again given by Eq. (90), and the 
last term in Eq. (90) vanishes in the limit N -^ oo (see 
the discussion at the end of Appendix C). 

Therefore, while not strictly stable, the system is 
asymptotically stable. Collecting results, we conclude 
that the second-order system is 



1 



-i»-. 


(94) 


wv + 1 1^1 W-, 

2 UJx^y^z 


(95) 



v- 



where the penalties represent boundary values, the 
penalty on the first equation is applied throughout the 
interior of the domain via Eq. (93), and the normal vector 
N is defined as in Eq. (76). 
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FIG. 3: Loo error of a plane wave ^{x,t) = sin(k-x — cut) 
evolved using the second-order system (solid green). Results 
of equivalent first-order evolutions are plotted for compari- 
son (dashed blue) . The domain consists of 27 identical cubic 
subdomains covering the region x,y,z £ [—15,15], and the 
successive resolutions have 5,7,9, and 11 Legendre-Lobatto 
grid points per subdomain along each dimension. In this test, 
k = (.3, .2, .1) and to — |k|. The Loo error is a moving average 
over an interval At = 50, which includes 50 data points. 



Numerical Tests 

The three-dimensional wave equation in Eqs. (94)-(95) 
with bulk penalty given by Eq. (93) is found to be robust, 
stable, and convergent in all of our tests. We have run 
simulations on multiple spherical shell, cylindrical shell, 
and cubic subdomains. As an example. Fig. 3 shows the 
Loo error (maximum nodal error) of a sinusoidal plane 
wave propagating through a domain consisting of 27 cu- 
bic subdomains. 

In this example, only the boundary-face part of the 
bulk penalty in Eq. (93) is used, and normal vectors are 
defined as in Eq. (76). Empirically, we find that the 
bulk penalties associated with edges and corners are not 
needed in this example. The incoming mode at an in- 
terface boundary is supplied by the adjacent subdomain, 
while at outer boundaries it is computed from the an- 
alytical solution. Time-stepping is performed using an 
explicit fourth-order Runge-Kutta method. The results 
of an equivalent first-order evolution are plotted for com- 
parison. 

A few empirical observations are worth noting. In 
practice, we find that the bulk penalty terms arising from 
edges and corners in Eq. (93) are not actually necessary 
for obtaining stable, convergent evolutions. In all the 
tests we have performed for scalar waves in flat space, 
the terms due to the faces of the boundary are sufficient. 
However, the additional terms in Eq. (93) may need to 
be included for complicated domain decompositions or in 
curved space applications. 

For more general systems of quasi-linear wave equa- 
tions (such as Einstein's equations in generalized har- 



monic form [7]), we find that it is sometimes necessary 
to include a boundary term enforcing continuity of the 
field ij} in the penalty. That is, one makes the replace- 
ment 6U- — > 6U- + Sip in the penalties. In the tests that 
we have performed, this is not required for a simple wave 
equation in fiat (or curved) space. 

An alternative to defining unique normal vectors on 
corners and edges is to use a so-called multi-penalty 
method. With a multi-penalty method, boundary con- 
ditions (and hence penalties) on edges and corners are 
defined to be the sum of those from the adjacent bound- 
ary faces. While this has the advantage of avoiding some 
of the issues with corners and edges, it makes obtaining 
analytical results such as Eq. (90) more difficult. Al- 
though we have not yet fully tested this alternative in 
curved space applications, we find that the multi-penalty 
method performs equally well for scalar waves in fiat 
space. 



IV. WAVE EQUATION ON CURVED 
BACKGROUND 



In this section we consider the application of the new 
penalty method to the evolution of a scalar wave on a 
fixed, curved background spacetime: 



V^V^tA = 0, 



(96) 



where V^ is the four-dimensional covariant derivative. In 
rewriting this equation as a first-order system, we use the 
standard 3-1-1 splitting of the metric: 



ds^ 



^a^dt^ + 7,y(dx' -I- P'dt){dx^ + (i^dt), (97) 



where a is the lapse function, /3* is the shift, and ^ij 
is the three-dimensional metric intrinsic to the constant 
time spatial hypersurfaces. It is assumed that a > and 
that the three-metric ^ij is positive definite. 

The wave equation in Eq. (96) can be rewritten in a 
standard way [20] as the first-order system 

i)= -aTi + I3'd^^p, (98) 

77= - aY^d^(l)j + I3'd^7r + aKn + aJV^, (99) 

0, = - ad.TT - TTd^a + ^kdrP^ + /3^9fc<^,. (100) 

Equation (98) is just the definition of the variable tt. As 
usual, the spatial derivative variable is defined as 



0j = diip. 



(101) 



The quantities K and J* in Eq. (99) are purely functions 
of the background spacetime: 



K 

J' 



1 



a-y 



1/2 



a7 






(102) 
(103) 
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where 7 = detjij. In deriving Eq. (100), the equivalence 
of interchanging indices in 



9i4'j = dj(t)i 



(104) 



has been assumed. This reduction to first order has there- 
fore introduced two constraints to the system: d = Cij = 
0, where 



Ci = 0j - d^ip, 
Cij = di(l)j - dj(j)i. 



(105) 
(106) 



The second-order in space equations are 

^ = - aTT + P^dii^, (107) 

TT= - aj'^d.dji; + f3'd,TT + aKir + aJ^diip. (108) 

This system avoids the introduction of the constraints 
in Eqs. (105)-(106) as well as the third set of evolution 
equations in Eq. (100). 

The characteristic variables and speeds of the second- 
order system are the same as those of the equivalent first- 
order reduction with (j)i — > diip: 

U^^^, Ao = -n'=/?fe (109) 

U±^T:±n'di^p, X±^±a-n'' 13k, (110) 

Ul> ^d,i> - n,n^dji>, Ao = -n'=/?fe, (111) 

where rf is the outward-directed unit normal vector to 
the boundary of the three-dimensional spatial domain. 
These are the same as the characteristic modes of the 
scalar wave in flat space in Eqs. (68)-(70), with mod- 
ified characteristic speeds. As discussed in Section III 
above, the "zero-speed" modes t/? can be considered to 
have arbitrary speeds in the second-order system [10]. 
The speeds Aq given above are chosen to be the same 
as those of the corresponding first-order system. Ad- 
ditionally, these are the coefficients that appear in the 
boundary fiux of the energy, and in this sense they are 
the preferred choice. 



Continuum Energy Estimate 

The energy density for this system is the same as for the 
flat space scalar wave in Eq.(72): 



1 



(tt^ + d'^d^ip). 



(112) 



The energy flux is found by computing the time deriva- 
tive of the energy 



In 



//'d'3 



(113) 



where fl is the spatial domain under consideration and 
7^'^d^x is the volume element. In addition to a bound- 
ary fiux, differentiating the energy gives rise to volume 



terms that depend on various derivatives of the back- 
ground {dia, diPj, or dijjk)- However, these volume 
terms can all be bounded by multiples of the energy it- 
self (or neglected entirely in the constant-coefficient ap- 
proximation), which is all that is required for proving 
well-posedness. One therefore obtains 

E<-1[ F'^a'^/^cfx + kE, (114) 

4 Jdn 

for some constant fc > 0. The flux integrand is 

F" ^\^Ul+ \+Ul + 2Ao L/°^C/]', (115) 

and the element of area in Eq. (114) is a^^^d^x, where 
a = det aij and (Jij is the intrinsic metric on the bound- 
ary surface. The continuum problem is therefore well- 
posed with boundary conditions that control incoming 
modes (those with A < 0). For a timelikc boundary, a 
boundary condition is needed on U- and possibly on C/f , 
depending on the sign of Ao. For a spacelike boundary, 
either all modes are incoming, or all modes are outgo- 
ing and no boundary conditions are required (e.g. on an 
excision boundary inside the horizon of a black hole). 

We could also have included a term a^tp'^ in the energy 
density, replacing Eq. (112) by 



This would give an additional term in E: 



(116) 



a^^pij-f^^^d^x^ / a^'iP{l3'd,'ip~aTr)j^/^d^x. (117) 
o Jn 

Integrating by parts in the first term on the right-hand 
side yields 

^ / n,p^^'a'/'d'x-^ f ^^d,{P^^^'^)d^x. (118) 
^ Jdci ^ Jn 

The latter term in this expression can be bounded by a 
multiple of the energy, while the first term contributes 
to the boundary flux. It may seem, then, that including 
the term a^i/'^ in ttis energy density would require the 
flux F" of Eq. (115) to be modifled. However, the entire 
right-hand side of Eq. (117) can in fact be bounded in the 
volume. Making use of the relation {/3^diip)'^ < d^ipdiip, 
we flnd 



a^ij{(3'd,i^-aTT)j^^^d^x 



< a(amax + |^|max)-B- 



(119) 



The addition of a term a^ip'^ to the energy density there- 
fore requires the constant k in Eq. (114) to be modi- 
fled, but not the flux F". Consequently, our conclusions 
about wcU-posedness and boundary conditions remain 
unchanged. It is interesting to note, however, that the 
same does not hold for the flrst-order system of Eqs. (98)- 
(100), because the first-order energy corresponding to 
Eq. (116) controls (j)i, but not ditp (and therefore the 
inequality in Eq. (119) does not follow). 
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Semi-discrete Energy Estimate 

The penalties in the semi-discrete equations need to be 
slightly modified from those of the flat space scalar wave 
system in Eqs. (94)- (95). To see how, consider the semi- 
discrete equations corresponding to Eqs. (107)-(108) with 
penalty functions p, q: 

ijj, ^ -aTTi + ...+p, (120) 

^^^^ -aJ'''^,^ki^^ + ... + q. (121) 

As usual, is it to be understood that the fields repre- 
sent grid values (e.g. ip = ''Pijk), and differentiation is 
implemented, for example, by matrix multiplication. For 
simplicity, we will assume that the physical domain un- 
der consideration has been mapped to the cube i7 with 
x,y,z G [—1,1]. We will also assume for now that the 
boundary is timelike, with C/_ the only incoming mode. 
As in the fiat space scalar wave system, the penalties 
will thus be proportional to 6U- = C/jJ?*^ — C/_. The 
semi-discrete energy is 



E=^[{7r,n) + {d'^,dr 



(122) 



where the discrete inner product is now defined by, for 
example 



(7r,7r) 



E2 1/2 
ijk 



TT^J^/^d^X. 



(123) 



Because of the presence of 7^'^ in the volume element, 
the quadrature integrals we encounter will in general no 
longer be exactly equal to the continuum integrals. 

The time derivative of the semi-discrete energy in 
Eq. (122) separates as usual into a continuum- like part 
plus a contribution from the penalties: 

^ -^continuum i -^penalties? V / 

where the penalty contribution is 

^penalties = (tT, q) + (9V, ^^p) . (125) 

Assuming the penalty q is defined as in Eq. (95) except 
for an overall factor go, we have 



(^,g) = i^|N|go7r,5C/_7i 



/2 



(126) 



ao 



where |N| is the magnitude (now with respect to jij) of 
the normal vector defined as in Eq. (76). The second 
term in Eq. (125) gives 

(9*V,5.p) = Y. \^\n'^^^Pp-f^/^ - (V,V^'^,p), (127) 

90 



where n* is the unit normal vector to the boundary, and 
Vj is the three-dimensional covariant derivative associ- 
ated with 7ij. With the penalty function p constructed 
in the volume according to Eq. (93), the last term in 
Eq. (127) asymptotically vanishes as in the fiat-space 
case, and we will therefore neglect it. Assuming that 
p has the same value as in Eq. (94) apart from an overall 
factor pq , it follows that 



(av, d,p) - -J 51 1^1 "'^^^ p^ ^^- ^' 



/2 



(128) 



dfl 



If we choose go = Po, then Eq. (125) for the penalty 
contribution to E becomes 



i^penalties = l^Yl '^'^O USU- 7^ 



/2 



(129) 



on 



Setting Pq ~ |A_|, we obtain the semi-discrete energy 
estimate 



i?<-2EiN|^"^'^' + ^^' 



(130) 



an 



for some constant fc > 0. The fiux integrand is 

F" = A_ (U^^^ - SUl) + X+Ul + 2XoU"'U°, (131) 

which resembles the continuum result of Eq. (115) with 
the addition of the negative term proportional to the mis- 
match of characteristic modes 6Ul. The sum over the 
boundary can be rewritten as 



E|n|7^^'(-) = E|n|e't^^'(-), 



(132) 



ao 



dn 



where |N|e is the magnitude (with respect to a Euclidean 
metric) of the normal one- form corresponding to N, and 
( • ) represents any integrand. In this form, the similarity 
to the surface integral 



{■)a'^^d^x 



(133) 



an 



in the continuum result of Eq. (114) is evident. 

We have assumed that C/_ is the only incoming mode, 
and in this case Eq. (130) shows that the semi-discrete 
system is asymptotically well-posed. If Aq < 0, then 
the boundary flux in Eq. (131) implies that a boundary 
condition is required to control U^ as well. Although we 
have been unable to see how to do this with penalties, we 
have found empirically that it is unnecessary to impose 
any boundary conditions on this mode; it is sufficient to 
enforce the condition on the incoming mode U- . 

For a spacelike boundary with all characteristic modes 
outgoing, no boundary conditions and hence no penal- 
ties are required. In that case, the boundary term in 
Eq. (130) is strictly non-positive. On the other hand, if 
the boundary is spacelike with all characteristic modes 



incoming, a boundary condition can be enforced on U- 
and C/+ by setting 



1 



INI 



1\ 



on 



5U-), (134) 

(|A+|(5[/+ + |A_|5L/_). (135) 



2 UJx UJy UJ. 

In this case, the flux integrand in Eq. (130) would be 



^ri 



+ 2XnU°'U?. 



sui 



(136) 



In summary, the second-order system with pcnahies is 



^ 



■SU-, 



7T = — aj^^didj^/j + . . . + 



|A_ 



INI 



UJx U->y UJz 



SU-, 



(137) 
(138) 



where as usual the penalties represent boundary values, 
the penalty on the first equation is applied throughout 
the interior of the domain via Eq. (93), and the normal 
vector N is defined as in Eq. (76). Furthermore, it is to 
be understood that the penalties with SU^ are applied 
only when [/_ is an incoming mode, and in the event that 
U+ is also incoming, the penalties are modified according 
to Eqs. (134)-(135). 



Numerical Tests 

In order to compute the error with respect to an an- 
alytical solution, we consider the inhomogeneous wave 
equation 



V. V^^ = 5, 



(139) 



where the source S is computed by substituting an an- 
alytical solution for tp into the left-hand side. With 
a source term, Eq. (107) remains unchanged, while 
Eq. (108) is modified by adding a term aS: 



',S. 



(140) 



As an example of a test problem, we use the following 
background metric (the Schwarzschild solution in Kerr- 
Schild coordinates): 



ds^ = - ( 1 + ) dt^ - 

AM , n 

H drdt + r^dil^ 

r 



2M\ , 
1 + ]dr^ 



(141) 



and consider an analytical solution of the form 

^'analytical = COs(c^ t)e~^^-'-«^" /^" Yi^. (142) 

We evolve the second-order equations Eqs. (137)-(138) 
with source term on a domain consisting of three concen- 
tric spherical shell subdomains. The innermost boundary 
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FIG. 4; Loo error of a scalar wave given by Eq. (142) on a 
Schwarzschild background in Kerr-Schild coordinates evolved 
using the second-order system (solid green). The results of 
equivalent first-order evolutions are plotted for comparison 
(dashed blue). The domain consists of three concentric spher- 
ical shells with radial boundaries at r = 1.9, 11.9, 21.9, 32.9 
(in units of M). The radial and angular resolutions {Nr,L) 
of the runs are (8,4), (14,6), (20,8), and (26,10). In this 
test, the following values for the analytical solution were used: 
ro = 17M, a = 2M, lv = 0.5M"\ / = 3, and m = 1. The 
Loo error is a moving average over an interval At = 25, which 
includes 50 data points. 



is placed just inside the horizon (located at r = 2M), so 
that no boundary condition is required there (all charac- 
teristic modes are outgoing). At interfaces between two 
subdomains, the value of the incoming mode U^'~^ is sup- 
plied by the adjacent subdomain, while on the outermost 
boundary, t/^^ is computed from the analytical solution 
in Eq. (142). Time-stepping is performed using an ex- 
plicit Runge-Kutta method. 

For a spherical shell subdomain, we employ a spec- 
tral basis composed of Legendre polynomials in the ra- 
dial direction scaled to the appropriate radial extent and 
spherical harmonics for the angular directions. The nu- 
merical approximant for particular truncations iV^ and L 
is therefore given by: 



N,. 



+1 



^ = 5ZE 5Z '^^lmP^{r)Yl,r,ie,(|)), 



(143) 



i=0 1=0 m=-l 



where Pi (r) represents the appropriately scaled Legendre 
polynomial, and aum are the spectral coefficients. Fig- 
ure 4 shows the Loo error in ip for this test problem as a 
function of time for several resolutions. For comparison, 
the results from an equivalent first-order evolution are 
plotted as well. 

The homogeneous second-order system of Eqs. (137)- 
(138) is also stable and convergent in tests that we have 
run with arbitrary initial data on a variety of back- 
grounds. For example, we have run simulations on 
a Schwarzschild background in Kerr-Schild, Painleve- 
GuUstrand [21], and fully harmonic coordinates, and on 



a Kerr background (with a spin up to a = 1) in Kerr- 
Schild coordinates. In those cases for which we do not 
have an analytical solution to supply a condition on U- 
at the outer boundary, we find good results by using a 
Sonimerfeld condition, assuming a solution of the form 



V' 



fit - r) 



(144) 



This translates into a condition on the incoming mode at 
the outer boundary: 



U 



BC 



v- 



(145) 



For the first-order system, as discussed below Eq. (38), 
one can use Chebyshev polynomials instead of Legendre 
for the radial basis in Eq. (143), and the results are com- 
parable (the same holds for Einstein's equations as well). 
In the second-order evolutions, while it is acceptable to 
use a Chebyshev basis for flat space applications, we find 
that the error is significantly larger (almost two orders of 
magnitude) in the curved background case. Note, how- 
ever, that we are not addressing the Chebyshev-Legendre 
method discussed below Eq. (38) here, but rather the use 
of a Chebyshev basis without modifying the index values 
of the penalty functions. Presumably, the Chebyshev- 
Legendre method would perform equally well, but we 
have not explored this modification. 



DISCUSSION 



13 



the sense that for a given accuracy goal, a smaller reso- 
lution is required in the second-order system. However, 
our focus has not been on comparing or analyzing code 
efficiencies, but rather on establishing the viability of a 
second order in space spectral method. We believe the 
second-order system has the potential to show substan- 
tial increases in efficiency for more complicated systems 
than those considered here. 

All of the energy arguments in this paper have assumed 
a grid structure that can be mapped to a cube. While 
it is straightforward to apply these methods to spherical 
or cylindrical shells, it is assumed that any dimension 
with boundaries (e.g. the radial dimension on a spherical 
shell) has a collocation grid that contains its endpoints 
(Gauss- Lobatto grid). For a domain containing the ori- 
gin, such as the unit disk, it is typical to use a radial grid 
of Gauss- Radau points, so that the endpoint at the origin 
is omitted (see e.g. Ref. [22]). We have not considered 
the generalization to such domains. 

We have shown how to evolve a multi-domain second 
order in space wave equation stably using spectral meth- 
ods. For more general systems, energy arguments like 
those given in this paper cannot be carried out. Never- 
theless, in the most important case we consider, namely 
Einstein's equations in generalized harmonic form, these 
methods work well. The reason is that the principal part 
of the equations is directly analogous to the scalar wave 
equation on a curved background [7]. We will report on 
these extensions in a subsequent paper [14]. 



We have not found a significant difference in effi- 
ciency between the first- and second-order forms of the 
equations for the simple systems considered in this pa- 
per. Even with equal time steps, the rates of the first- 
and second-order codes (using explicit time-stepping) are 
comparable (within about 10% of each other), despite 
the absence of the third set of evolution equations (the 
0i equations) in the second-order system. One reason for 
this is that second derivatives are computationally more 
expensive than first derivatives on arbitrary domains, be- 
cause of the coordinate transformations involved. Sym- 
bolically, the transformation of a first derivative to new 
(barred) coordinates involves a multiplication by the Ja- 
cobian J of the transformation: 



dip ~ J dip, 



(146) 



whereas the transformation of a second derivative re- 
quires the Hessian H and the first derivative as well: 



d'^ip = JJd^ip + Hdip. 



(147) 



It is evident that the second-order evolutions generally 
have smaller errors than their first-order counterparts at 
given resolutions by as much as two orders of magnitude, 
as can be seen in Figs. 3-4. To a degree, then, this makes 
the second-order evolutions somewhat more efficient in 
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Appendix A: Gauss-Legendre-Lobatto Quadrature 

Here we provide some of the properties of Gauss- 
Legendrc-Lobatto quadrature [6]. The basis functions 
on [— 1 , 1] are the Legendre polynomials P„ (x) . This is 
a convenient choice for obtaining analytical results be- 
cause the Legendre polynomials are orthogonal under a 
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weighting function of unity: 



p{x)Pn{x) P,n{x) dx = (5„„, 



with p(x) = 1. The (N + l)-point quadrature rule, 



N 



(Al) 



/ u{x) dx ~ 2. ^i '^'{Xi) 



(A2) 



is exact if u{x) is a polynomial of degree 2N — 1 or less. 
The 7V+1 nodes Xi are 



Xq -- 


= -1, 


(A3) 


XN -- 


= + 1, 


(A4) 


Xi -- 


= the roots ofP^(a:;) 


(A5) 




for < 2 < iV, 


(A6) 



and the weights coi are given by 

2 



N{N + 1)[Pn{x,)Y^ 



(A7) 



Note that there is no known explicit formula for the roots 
of P'j^{x) — they must be found numerically. A function 
i/'(a;) is approximated by an 7V*^-order interpolating poly- 
nomial 'tpN{x), which can be expressed as 



N 



iPn{x) = '^■^p{xi)C^{x), 



(A8) 



i=0 



where Ci{x) are cardinal functions satisfying Ci{xj) 



They can be written as 



Q{x) 



-{l-x^)P'^{x) 



N{N+1)Pn{x,){x~x,)' 



(A9) 



Differentiation can be computed via matrix multiplica- 
tion from 



N 






(AlO) 



i(i) 



where D)',' = CAxA is the first-order differentiation ma- 
trix. The second-derivative matrix is defined similarly 
and satisfies I?^^^ = D^^^D^^'. An efficient algorithm 
for computing pseudo-spectral differentiation matrices is 
given in Rcf. [5]. 

If /, g are two iV^'^-order polynomials, summation by 
parts follows naturally because the product fg' is a poly- 
nomial of order 2iV — 1 or less: 



N 



(/,5')=E^^/^5: = /^ 



U=0 



(/',3)- (All) 



Summation by parts generalizes to higher dimensional 
inner products in a straightforward way. For example, if 
/ and g are 2-d polynomials in x and y: 



N 



i,j=0 

N 

= ^'^]{f9)\,=a- {.f.9^9)■ 



iA12) 
(A13) 



3=0 



Appendix B: Proof of Inability to Generalize ID 
Penalty Function 

In this section we will show that in two or more dimen- 
sions the inner product 



{d'd,i^,p) 



(Bl) 



that arises in the energy arguments discussed in this pa- 
per cannot be made to vanish in general with a penalty 
function p that satisfies the boundary conditions. We 
will argue by counting degrees of freedom. For simplic- 
ity, consider the two-dimensional case and let the domain 
be a square with A^-|-l grid points along each dimension. 
Instead of using a basis of Legendre polynomials, con- 
sider a (non-orthogonal) basis of functions x*y^. A 
scalar field is thus approximated on the grid as a two- 
dimensional interpolating polynomial of the form 



■0 = ^ aijx'y^. 

a<i,j<N 



(B2) 



There are (A^ + 1)^ basis functions and hence the same 
number of degrees of freedom in the function -0. The 
penalty function p must satisfy 4A^ boundary conditions 
on the square. 

Now consider operating on the expansion of V in 
Eq. (B2) with the Laplacian d^ + 9y. The effect of this 
operation on a term x'^y^ is essentially 

x'y^ ^x'-^y^ +x'y^-^. (B3) 

Since we are only interested in counting the degrees of 
freedom that remain in V^i/", we only need to retain one 
of the terms in Eq. (B3): 



x^'y^ -^ x' ^y\ 



(B4) 



j=0 



By doing this, we will at worst undercount the degrees 
of freedom in V^ V- This leaves terms of the form x*~^j/^ 
for 2 < i < N and < j < N, which implies that there 
are at least (A^ + 1)(A^— 1) degrees of freedom remaining 
in the Laplacian. 

There are thus at most (N+l)^ - {N+1){N-1) = 2N+2 
degrees of freedom for constructing a penalty function 
that is orthogonal to V'^ip (for arbitrary -0), which is not 
enough to satisfy the 4A^ boundary conditions. The same 
argument can be applied in any number of dimensions. 
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In particular, in the three-dimensional case we find that 
there are at most 2{N + 1)'^ degrees of freedom for con- 
structing the penalty function — not enough to satisfy the 
6iV^-|-2 boundary conditions, which proves the assertion 
made below Eq. (90). 



Appendix C: Derivation of 3D Penalty 

In this section the form of the three-dimensional bulk 
penalty given by Eq. (93) will be derived. The goal is to 
minimize the inner product 



{d'd.i'^p), 



(CI) 



with the values of the penalty function p on the boundary 
given. First, we will revisit the one-dimensional problem 
on the interval [—1,1] from a new point of view. In Sec- 
tion II C it was shown that the one-dimensional inner 
product {ip" ,p) vanishes when p is constructed from the 
functions f,g defined in Eqs. (59)-(60). Recall that the 
functions / and g were constructed from Pjy and Pn-i 
to be orthogonal to ip" . 

Let us start over and consider the one-dimensional 
penalty function p to be unspecified, except on the 
boundaries. Suppose also that there is no boundary 
condition at xjv — +1, so the penalty function satisfies 
Pn — 0. The boundary condition at xq = —1 fixes the 
value po, and we can view the values pi for < « < iV as 
free parameters for minimizing the inner product: 



N-l 



{ip",p) = cjo-^o Po + X! '^^V'j'k- 



(C2) 



The case with a boundary condition at x = +1 and 
not at X = —1 (so Pq = 0) is similar, and we conclude 
that the interpolation weights CAr(xi) for approximating 
a function at xn = +1 based on its values at points Xi 
for < i < N are given by 



CN{xi) = -gi 

ujn 



(C6) 



where g, are the grid values for < i < N oi the function 
g defined in Eq. (60). 

Now let us consider the two-dimensional problem on 
the square [—1, 1] x [—1, 1]. In the following we will use 
the index i exclusively for summing over x values and 
j for y. Our goal is to construct the values of p on the 
interior of the domain so as to minimize the inner product 



{Axp,p) = y^^uj^ujjAxpijp^j, 



(C7) 



where A represents the Laplacian operator d^di, and we 
consider the values of p on the boundary to be given. 
Consider a point on the edge at (xo,yj), for example. 
The term in the inner product due to this point is 

ojQUjjAipojPoj. (C8) 

Now define p on the interior along the j*'* row to be 

Pij^ Poj.fi, (C9) 

just as in the one-dimensional case. Using the identifi- 
cation of / as interpolation weights from Eq. (C5), the 
contribution to the inner product from the interior of this 



N-l 



{^'>P,P)\jt^row = Yl ^^^j^'^ijP^3 



(CIO) 



This will vanish if and only if 

N-l 



N~l 



c = E 






(C3) 



where it is safe to assume po ^ (if po = 0, then p = 
as there would be no need for a penalty function). Since 
ip'^x) is an arbitrary {N—2)-OTdcT polynomial, this equa- 
tion defines the ideal interpolation weights co(xi) for ap- 
proximating a function at xq = —1 based on its values 
over a stencil of points x, for < i < A'^: 



Af-l 



< = J2 ^0{x^)i^^ 



(C4) 



i=l 



Assuming the grid points are Gauss-Legendre-Lobatto 
points, we can therefore make the identification 



co[xi) = = }i, 

ujqPo ljo 



(C5) 



where fi are the grid values for < i < iV of the function 
/ defined in Eq. (59), and we have used the fact that 
Eq. (C3) holds when the penalty p is defined by Eq. (61). 



E^^^. 



i=l 



(Cll) 



~ - wqWj AV'ojPoj, (C12) 

which approximately cancels the term from the point on 
the edge in Eq. (C8). In Eq. (Cll) we have written coi 
for the interpolation weights co(xi) defined in Eq. (C5). 
Next, consider a point at a corner, say {xN,yo)- The 
term in the inner product due to this point is 

WATWoAV'TvoPwo- (C13) 

Define p on the interior of the domain to be 

Pi] = -PNo9i.fj- (C14) 

The contribution of pij to the inner product on the inte- 
rior of the square is now 



{Atp,p)\ 



N-l 

E 



uj^ujjAipijPij 



(C15) 



■ ujn(^oPno 



N-l 

E 



Aipij CNiCoj (C16) 

■UJNi^oAlpNOPNO, (C17) 
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which approximately cancels the contribution from the 
point on the corner in Eq. (C13). Following this proce- 
dure, we construct p on the interior by adding a contribu- 
tion from each boundary segment: 4 edges and 4 corners 
on this 2-d domain. Explicitly, this gives 



Pij ^Pojfi + PNj gi + Piofj + PiN gj 
- Paofif] - PoNfigj - Pno gifj 

-PNNgigj- 



(C18) 



This generalizes to three or more dimensions in a 
straightforward way. Each term in pij has a number of 
products of / or (7 equal to the codimension of the bound- 
ary segment it depends on. The only caveat is that the 
sign of the terms added to p should be (—1)™+^, where m 
is the codimension of the boundary piece producing the 
term. This is evident in the 2-d example above where the 



terms in Eq. (C18) due to the corners are negative. The 
sign difference arises simply because of the negative sign 
in the relation between the interpolation weights cq,cj\j 
and the functions f,g in Eqs. (C5)-(C6). 

With the penalty function p constructed according to 
the above procedure, the inner product of p with any 
analytic function h (hence Aip) satisfies 



{h,p) -> 0, as N —)■ oo. 



(C19) 



In particular, we have shown that the last term in 
Eq. (90) asymptotically vanishes, as claimed below 
Eq. (93). Moreover, while we have not bounded the error 
for a given resolution, the inner product in Eq. (C19) will 
be as small as possible in the sense that it vanishes for 
the polynomial approximations to h up to order N — 2. 
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